Universal features and tail analysis of the order-parameter 
distribution of the two-dimensional Ising model: An entropic 

sampling Monte Carlo study 

Anastasios MalakiJ*] and Nikolaos G. Fytas 
Department of Physics, Section of Solid State Physics, 
University of Athens, Panepistimiopolis, 
GR 15784 Zografos, Athens, Greece 
(Dated: September 15, 2008) 

Abstract 

We present a numerical study of the order-parameter probability density function (PDF) of the 
square Ising model for lattices with linear sizes L = 80 — 140. A recent efficient entropic sampling 
scheme, combining the Wang-Landau and broad histogram methods and based on the high-levels 
of the Wang-Landau process in dominant energy subspaces is employed. We find that for large 
lattices there exists a stable window of the scaled order-parameter in which the full ansatz including 
the pre-exponential factor for the tail regime of the universal PDF is well obeyed. This window 
is used to estimate the equation of state exponent and to observe the behavior of the universal 
constants implicit in the functional form of the universal PDF. The probability densities are used 
to estimate the universal Privman-Fisher coefficient and to investigate whether one could obtain 
reliable estimates of the universal constants controlling the asymptotic behavior of the tail regime. 
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I. INTRODUCTION 



A significant achievement in the theory of equilibrium critical phenomena was the 
confirmation of universality and scaling hypotheses and the calculation of critical expo- 
nents The finite-size scaling technique has proven to be an extremely 
reliable and powerful method for determining the critical properties of low-dimensional sys- 
tems, and several review articles have appeared covering the original finite-size scaling theory 
and later advancements [?, [s, y, [l^. In practice all universality claims have been put to 
several tests using data obtained from numerical simulations of finite systems and the extrac- 
tion of universal behavior from such studies remains today of vigorous research interest. Of 
particular interest is the view that, a strong hallmark of a universality class may be obtained 
via the probability density functions of the main thermodynamic variables of the system at 
criticality llj. In finite systems the order-parameter may be characterized by a probability 
distribution function (PDF) and this distribution may be scaled appropriately to provide a 
sequence of densities that appears to converge rapidly to a unique scaling PDF 12 , [isl, [l^ . 
We may think of this limiting scaling function as the key ingredient specific to a universality 
class, like the critical exponents 15| . This idea has been elaborated in recent years by many 
authors, several papers are devoted to the estimation of such universal distributions, and 



the subject is of great and growing interest 
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Most properties of finite-size scaling functions are known from numerical simulation of 



critical systems 12 




from field theoretic renormalization group calculations 



241 Analytical results originate 



26|, |27|], conforma! 



but also from a generalized classification theory of phase transitions 
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ield theory 



3, 



31 



In 



a few cases some of the analytical predictions seem to have been confirmed by numerical 



simulations |. A detailed investigation of the tail regime was carried out in one of 

the early multicanonical Monte Carlo simulations by Smith and Bruce {2^ for the two- 
dimensional Ising model, using square lattices of size L = 32 and L = 64. Even though 
these authors managed to measure extremely small tail probabilities with high accuracy, 
and found signs for the far tail regime conjecture, the theme was not fully elucidated. It 
is well-known that the process of recording the very small probabilities in the tails of the 
distribution requires specialized numerical techniques and sufficient statistical accumulation 
is necessary in order to probe the tails confidently. 
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Some studies using traditional Monte Carlo simulations have attempted to estimate the 
universal density over the last years, but failed in establishing the true behavior at the far 
tail regime of the critical order-parameter distribution. Traditional Monte Carlo sampling 
methods have increased dramatically our understanding of the behavior of the standard 
classical statistical mechanics systems. The Metropolis method and its variants were, for 
many years, the main tools in condensed matter physics, particularly for the study of crit- 



ical phenomena 
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37l |. However, this standard approach has certain crucial 



weaknesses. Importance sampling is trapped for very long times in valleys of rough free 
energy landscape, as in complex systems in which effective potentials have a complicated 
rugged landscape. For large systems these trapping effects become more pronounced and 
the Metropolis method, but also its variants, become inefficient. Moreover, importance sam- 
pling methods are unsuccessful in recording of the very small probabilities in the tails of the 
critical order-parameter distribution. As was pointed out by Hilfer et al. 2l(] the study of 
the tail regime requires special numerical techniques 19, 20, SsJ. 

Several new efficient methods that directly calculate the density of states (DOS) of clas- 
sical statistical models play a dominant role in overcoming the above mentioned problems 



36, 



39|, 



in recent years. A few remarkable examples of such methods are the entropic 
multicanonical [40|, broad histogram (BH) [41I], transition matrix and Wang-Landau 
(WL) 43| methods. It is now possible to study more effectively the finite-size scaling prop- 
erties of statistical models by using entropic Monte Carlo techniques. The multicanonical 
Monte Carlo method has already been applied for the evaluation of the tail regime of the 
universal PDF of the order-parameter 211 ^^"^ present paper follows these attempts 
by using a new simple and efficient entropic technique. This technique will be referred to 



as the critical minimum energy subspace CrMES-WL entropic sampling scheme [4J] and is 
in fact a program for the simultaneous estimation of all thermal and magnetic finite-size 
anomalies of the statistical systern^The method was recently presented and successfully 
tested on the square Ising model [44]. The central idea is to optimize the simulations by 
using in a systematic way only the dominant energy subspaces appropriate to the finite 



system at the temperature range of interest. As shown in Ref. [4J] this restriction speeds up 
our simulations, but also gives a new route for critical exponent estimation by studying the 
finite-size scaling of the extensions of the dominant subspaces. This scheme is expected to be 
much more efficient than the Metropolis algorithm, as already pointed out in a comparative 
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study [4^, illustrating its superiority in the recordings of the far tail regime. Overall, the 
aim of this contribution is twofold: firstly to show that the CrMES scheme is sufficient for 
the study of the tail regime and therefore to propose a sensible optimization route. Sec- 
ondly, to yield new evidence confirming the tail regime conjecture and to provide an early 
estimation of the universal parameters involved. 

The rest of the paper is organized as follows. In the next section, we briefiy review the 
CrMES- WL entropic scheme, including the N-fold way implementation of this scheme. This 
approach is more efficient and makes available an additional approximation for the DOS, 
based on the BH method. In Sec. [TTllwe study the order-parameter distribution of the square 
Ising model and we verify its asymptotic behavior in the far tail regime. The relevance of 
our findings to the universal Privman-Fisher coefficient is also discussed and used to observe 
the self-consistency of our estimations. Our conclusions are summarized in Sec. IIVI 

II. ENTROPIC SAMPLING IN DOMINANT ENERGY SUBSPACES 

Let us briefly describe the Monte Carlo approach implemented in order to generate the 
numerical data used in the next section. Since the main ingredient of this approach, decisive 
for its efficiency, is the restriction of the energy spectrum called CrMES restriction, we 
shall also explain, although briefly, some key ideas about its implementation. For a full 
description, alternative deflnitions and further technical details of the CrMES scheme one 



should consult the original papers, where the method has been recently tested [4^, |45|, |46 |. 
The main idea is to produce accurate estimates for all flnite-size thermal and magnetic 
anomalies of a statistical system by using a DOS method, based on WL random walks, in 



an appropriately restricted energy subspace {Ei, E2). In particular it was shown 4J] that it is 
quite accurate to implement this restricted scheme and at the same time accumulate data for 
the two-parameter energy, magnetization {E, M) histograms, using for their recording only 
the high-levels of the WL diffusion process. At the end of the process the flnal accurate WL 
or the accumulated BH approximation of the DOS and the cumulative {E, M) histograms, 
are used to determine all properties of the statistical system. This approximation is accurate 
in a wide temperature range, depending on the extension of the dominant energy subspace 
used, around the critical temperature of the system and this range can be chosen in such a 
way that includes all pseudocritical temperatures of the flnite system. 
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A multi-range WL algorithm 43|] is implemented to obtain the DOS and the {E, M) 
histograms in the energy subspace {Ei,E2). The WL modification factor (fj) is reduced 



at the jth iteration according to: /i 



Jj-l, J 



2, J fin- 



The density of 



states obtained throughout the WL iteration process may be denoted by Gwl{E). This 
density and the high-level (j > 12, see also the discussion below) WL {E, M) histograms 
(denoted by Hwl{E, M)) are used to estimate the magnetic properties in a temperature 
range, which is covered, by the restricted energy subspace {Ei,E2). The N-fold version 



of the WL algorithm [47|, |48| makes available an additional approximation for the DOS 
(denoted by Gbh{E)), based on the BH method. We shall use only the approximation of 
Gbh{E) corresponding to the minimum energy transitions {/\E = 4) in the N-fold process. 
Gbh{E) is obtained from the well-known BH equation [41]: G{E){N{E, E+AE))e = G{E+ 
AE){N{E + AE,E)) E+A.E, where N{E, E + AE) is the number of possible spin flip moves 
from a microstate of energy E to a. microstate with energy E + AE. The micro canonical 
average of these numbers is estimated at the end of the CrMES-WL entropic scheme, with 
the help of the corresponding appropriate histograms recorded in the same high-level WL 
iteration range, used also for the recording of the {E, M) histograms. 

The probability density of the order-parameter at some temperature of interest T may 
be expressed as follows |4j] 

T.E^iE,,E,)[HwL{E,M)/HwLmGwL{E)e-P'' 



Pt{M) 



E 



E£(Ei,E2) 



{E)e- 



■(3E 



where 



Hwl{E) = J2Hwl{E,M) 



(la) 



(lb) 



M 



and the summation in M runs over all values generated in the restricted energy subspace 
{Ei,E2). Since the detailed balance condition depends on the control parameter it is 
suggested 4J| that only the high-level recipes or their N-fold versions should be used for 
recording the [E, M) histograms. This practice yields excellent estimates for all magnetic 
properties as has been shown in detail in Ref. 44 1. 

The performance limitations of entropic methods, such as the WL random walk and 
the reduction of their statistical fluctuations have recently attracted considerable inter- 
est |5ll, l52, l53|]. For the CrMES entropic scheme, presented here in Eg . Q , an extensive 

44| . In particular. 



comparative study using various implementations was presented in Ref. 
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this study has clarified the effect of the used range of the WL iteration levels on the mag- 
netic properties of the system and also the effect of one of the simplest refinements of the 
WL algorithm. Statistical fluctuations are reduced, as usually, by multiple measurements 



but also by using the separation refinement proposed by Zhou and Bhatt [52j, in which the 
WL DOS modification is applied after a number S {S = 16) spin-fiips. For the recordings 
of the {E, M) histograms, the range j = 12 — 24, of the WL iteration levels, was shown to 
give accurate estimates for moderate lattice sizes (L = 10 — 120), while a range of the order 
j = 18 — 28 would be a more safe choice for larger lattices. Furthermore, as it has been 



shown recently 4J, |52|, |53l] , the histogram fiatness criterion of the WL scheme for reducing 



the modification factor should be treated with caution, and in order to avoid strong statisti- 
cal fiuctuations in the final approximate DOS, enough statistics should be obtained in each 
WL iteration. 

Let us provide here some further technical details on the implementation of our entropic 
scheme. For a large system, we divide the total CrMES in several energy ranges of the 
order of 50 — 60 energy levels each, overlapping at their ends in at least 3 energy levels. 
When combining these energy intervals, simple averaging is applied on the entropies of two 
neighboring pieces at their overlapping end. Then, the DOS for each piece is adjusted to 
the DOS of its previous neighbor and continuing this procedure we finally obtain the DOS 
in the total CrMES. We follow a mixed WL process in which in the first stage (j = 1 — 11 or 
j = 1 — 17) we use the WL algorithm and in the second stage {j = 12—24 or j = 18—28) its N- 
fold version. This mixed process was also used in a previous study [48] , where a generalization 



oft 
al. 



le N-fold version was presented. In the present application only the original Schulz et 



471] N-fold version has been used (case c = 1 in Ref. 48] ). The accumulation of the 



microcanonical estimators, necessary for the application of the BH equation [41|, takes place 
only in the N-fold stage of the process. At these high levels of the WL random walk, the 
incomplete detailed balanced condition has not a significant effect on the BH microcanonical 
estimators and the particular multi-range approach seems to be optimal in both time and 
accuracy requirements. Besides the worse time requirements, our tests indicated also more 
significant effects (presumably coming from incomplete detailed balanced condition) for a 
multi-range approach using much larger energy intervals. Furthermore, an accurate run 
requires enough statistics to be obtained during each step of the process. 

The resulting WL and BH DOS approximations were tested using our accurate finite-size 



6 



data for the critical specific lieat and the asymptotic formulae discussed in Refs. |45l . |50 |. 
Using several independent random walks, we have verified that both the final WL and BH 
estimates for the critical specific heats (C(Tc)) are very accurate, with errors of less than 
1%, even for large lattices of the order of L = 200. In many cases, the BH estimates 
showed smaller statistical fluctuations when different runs were compared. However, in 
cases in which insufficient simulations were performed for large lattices {L = 140 — 200), 
a significant underestimation was observed for C(Tc) and this distortion was then much 
stronger for the BH estimates. A further interesting test, substantiating the statistical 
reliability of our approach, was performed for the lattice L = 50 by using the exact DOS 
obtained by the algorithm provided by Beale 54|. In this test, we applied a one-range exact 
entropic N-fold sampling to obtain an approximate BH DOS and this was then used to obtain 
an estimate C{Tc). When this one-range exact entropic estimate of C{Tc) was compared 
with the corresponding estimate obtained by our multi-range and approximate entropic 
scheme, we found excellent coincidence and in effect the same order of relative errors (0.001 — 
0.0001) compared to the exact result. Similarly, such small relative deviations, between 
the exact and the approximate entropic schemes, were also observed in the corresponding 
magnetic properties (for instance for a/ (m^)) and in fact these deviations were less significant 
than the fluctuations coming from limited magnetic sampling, as observed in different runs. 



Therefore, we suggest that the unified multi-range implementation of the Wang-Landau [43 1 
algorithm and the BH method of Oliveira et al. 4l| may be advantageous in studies using 
entropic schemes. 

In the original paper 45|, the CrMES method was presented by restricting the energy 
spectrum around the value producing the maximum term in the partition function at 
some temperature of interest. The restriction was imposed by requesting a specified relative 
accuracy (r) on the value of the specific heat. In other words, the restriction of the energy 
spectrum produces an error on the value of the thermodynamic parameter, at the particular 
temperature, and r measures the produced relative error. This relative error is set equal to a 
small number (r = 10^^), which is, in any case, much smaller than the statistical errors of the 
Monte Carlo method (the DOS method) used for the determination of the thermodynamic 
parameter (for instance the specific heat). It was also shown that a systematic study of the 
finite-size extensions of the resulting dominant subspaces produces the thermal exponent 
a/u with very good accuracy As pointed out above, this idea may be simultaneously 
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applied [44] for all finite-size anomalies including the magnetic anomalies determined from 
the {E,M) histograms. Furthermore, very good estimates of the critical exponent j/u 
have been obtained [44!] by studying the finite-size extensions of analogous critical minimum 
magnetic subspaces (CrMMS) defined also below. In this case, the restriction on the order- 
parameter space around the value M that maximizes the order-parameter density at the 
critical temperature Tc is imposed on the probability density function, as shown in Eq. ([2]). 
That is, the location of the dominant magnetic subspaces (M_,M_|_) may be obtained by 
comparing the end-point densities with the peak-height of the distribution 

~ M^<, (2) 
PtJM) 

and the critical exponent 'j/u may be estimated from the following scaling law (AM = 

(M+ - M_)) y _ 

(3) 

One should note here that, the finite-size extensions of the above defined CrMMS can be 
calculated by any Monte Carlo method producing the order-parameter distribution. One 
could as well implement the Metropolis algorithm to find estimates of the extensions involved 
in Eq. ([3]). However, this will yield a marked underestimation of 7/1/ as a result of the 
statistical insufficiency of this traditional algorithm in the tail regime. This effect was shown 
to be a result of the very slow equilibration process of the algorithm in the far tail regime of 
the order-parameter distribution 4J]. On the other hand, the CrMES-WL entropic scheme 
was shown to produce very good estimates of the critical exponent 7/z/ 4J] and this can be 
taken as an indication in favor of the suitability of this method in studies of the universal 
distributions and in particular for their tail regime. 

Let us now address the question of adequacy of the CrMES restriction to deal with 
the far tail regime of the order parameter PDF. In previous Monte Carlo studies 20, 121 1 
extensive simulations were undertaken in order to reach, as close as possible, the saturation 
regime {M/N ~ 1, = L^) of the order-parameter. This practice is related to the fact 
that the part of the order-parameter spectrum determining the tail regime is not formally 
known. Scaling the PDF introduces a variable x = m/ \/JrnF) and the (right) tail regime is 



expected to be detected in the range a; > 1 



13 



201 ] . Consequently, a main question concerns 



the sufficiency of the CrMES scheme. Does such a restriction on the energy space yield a 
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reliable approximation in the tail regime? In particular, we would like to know whether 
this restriction would allow us to simulate the order-parameter in the appropriate range to 
confirm the tail regime conjecture. The following observations provide strong evidence to 
this important question and our findings in the next section establish explicitly the fact that 
the implementation of a CrMES restriction permits an asymptotic evaluation of properties 
related to the tail regime. 

In Fig. [1] we illustrate the effect of various restrictions on the universal PDF (for details 
see next section) for a large lattice size L = 140. The simulation was carried out in a wide 
energy range {Ri in Fig. [1]) including the CrMES at the exact critical temperature {R2 in 
Fig.[T]). Specifically, counting the energy levels from the ground state with an integer variable 
ie{= 1,2,3, ...), where ie = 1 corresponds to the ground state and the energy of a level is 
E = —2N + 4- (ie — 1), the simulation was carried in the wide subspace Ri: ie = 1950 — 3800. 
The CrMES restriction applied to produce a relative accuracy (r = 10~^) on the critical 
specific heat (see the original paper on CrMES yields the subspace R2: ie = 2228—3514, 
while the restriction defined in Eq. ([2]) gives the subspace ie = 2246 — 3491, which by 
definition determines the extent we probe the tail of the order-parameter distribution. We 
note that this later subspace is included in the CrMES determined from the specific heat 
condition (using the same level of accuracy r) and that this is a general property which 
does not seem to depend on the lattice size, or even on the model as far as we had the 
opportunity to observe. This explains the striking observation from Fig. [T] that the PDF's 
obtained from the subspaces -Ri and R2 completely coincide in the x-range common in both 
subspaces. The difference between the two curves is always smaller than the accuracy level 
r and, as one can observe more clearly from the inset, the part of the i?i-PDF not accessible 
by using the i?2(CrMES) subspace is the range: x = 1.47 — 1.59. However, the i?i-PDF 
appears to be already flat in this range x = 1.47— 1.59 and fitting attempts will be unstable, 
producing large errors. Thus the saturation range is sensibly excluded from our simulations 
by the CrMES restriction. The third curve, denoted as R3 in Fig. [1] corresponds to a PDF 
obtained in the subspace R3: ie = 2620 — 3800, restricted severely from the saturation side. 
This curve is presented in order to observe that errors will be introduced if a restriction 
on the energy space is inappropriately applied. Noting that the errors in determining the 
end-points of the CrMES are of the order of 2 — 10 energy levels (at these lattice sizes) we 
conclude that critical minimum energy subspaces, with a reasonable accuracy level of the 
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order r = 10 ^ — 10 ^, will be sufficient for the study of the tail regime. 

III. UNIVERSAL FEATURES AND THE TAIL REGIME OF THE ORDER- 
PARAMETER DISTRIBUTION 

The universal scaling distribution of the order-parameter may be obtained from the mag- 
netization distributions Pm{iTL), constructed with the help of the numerical scheme outlined 



in the previous section, as follows 



20| 



p{x)dx ~ Pm{m)dm, x = m /v/(m2), m = M/N. (4) 

Figure [2] shows these distributions for lattice sizes L = 60 and L = 120. The data used were 
generated by the CrMES-WL(N-fold: 14-24) entropic scheme using a separation refinement 
S = 16, as explained in the previous section. The curves shown in this figure illustrate the 
densities obtained only via the BH approximation for the DOS and not the corresponding 
WL approximation of Eq. (jl]). This practice will be followed bellow in all our figures, except 
for the cases, indicated in the figures, where both the WL and BH DOS are used to construct 
and illustrate the corresponding approximations for the probability densities. Let us point 
out that the curves shown are lines that pass through all the points representing the sampled 
values of the order-parameter. For a large lattice, there will be several thousands of such 
points, since each of them corresponds to a possible value of m (= M/N, M = 0, 2, 4, N) 
of the finite system. The density of points on the a;-axis (m-axis) grows with the lattice 
size, and we should expect to having N/2 = L^/2 points in the positive x-axis (m-axis), 
provided that all energies and all corresponding order-parameter values were sampled by the 
WL process. However, as discussed in Sec. HIl the points corresponding to the saturation 
regime are not sampled, since this regime is excluded by the CrMES restriction applied on 
the energy spectrum. In order to illustrate the density functions in Fig. [2] we have chosen 
to identify their peak-heights to the same value, set in Fig. [2] equal to unity {p{x*) = 1, 
where x* the most probable value). This is equivalent to multiplying the universal PDF 
by a factor, which could in principle be weakly L-dependent and this dependence will be 
discussed further bellow. 

Let us now consider the main conjecture for the large-x behavior of the universal function 
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p{x) [19|, [2 

p{x) ~ poox'f' exp(-aooa;^^^), (5a) 

with 

^ = ^ (5b) 
and Poo, floo universal constants. The structure of the exponential ( |5all has been suggested by 



rigorous results for the two-dimensional Ising model 49] and is also consistent with Monte 



Carlo studies of the Ising universality class 



13| . The studies of Refs. [l9|, l20|] have provided 



some evidence for this conjecture and in particular for the prefactor and the relation of the 
exponent ip to the critical exponent 6. For the 2D Ising model the exponent should have the 



value ifj = 7, if, of course, the prefactor hypothesis is valid. Smith and Bruce [20| provided 
numerical support for this value, but their study was not completely conclusive since it was 
carried out only for relatively small lattices {L = 32 and L = 64) and the x-window in which 
the value ip = 7 was observed was actually quite narrow. We now present results for several 
lattice sizes (L = 80, 100, 120, and 140) reinforcing this conjecture in a very wide x-window. 
Following Smith and Bruce [20j we fix the exponent 6 in the exponential factor of Eq. and 
fit our results [x > 1) in x-windows, each one corresponding to 50 different magnetization 
values, sampled during the WL(N-fold: 12-24) process. Fig. [3] shows a very clear signature of 
the prefactor law (]5bp which upholds in a large x- window only for the large lattices {L = 120 
and L = 140 are shown). On the other hand, for smaller lattice sizes L < 80 (L = 80 is 

n 

shown) the picture is similar to that presented in Ref. [20[ and the expected value is obtained 
only in a small x- window. In Fig. Owe have illustrated the behavior of the estimates for the 
exponent ip, and for the large lattices {L = 120 and L = 140) both the WL and BH cases 
are shown. From this figure we observe that the window x = 1.2 — 1.3 is very stable for the 
sizes L = 120 and L = 140, and the estimates for the exponent ip remain close to the value 
ip = 7, beyond small fluctuations, for even larger values of x. It appears that this stable 
x-window may be very convenient for further fitting attempts. 

Encouraged from the above finding, we now use this stable x-window to perform further 
fitting attempts, allowing this time both exponents in the exponential and in the prefactor 
term vary, and use tp as a. free parameter by assuming the validity of Eq. fISbl) . In this way 
we examine whether the particular x-window provides a good independent estimation of the 
exponent 6, based on the full tail regime conjecture [Eq. ([5])]. In other words, ip, and 
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Poo are the three free fitting parameters in applying Eq. ([5]) to our numerical data. Figure H] 
depicts the three-parameter fitting attempts for L = 100 and L = 120 using both WL and 
BH approximations of the universal PDF. The values of the estimates for the exponent ip and 
their fitting errors are illustrated in this figure. The estimates obtained for the exponent ijj 
are quite good in all cases and their average = 7.041 is accurate to two significant figures, 
which is a rather pleasing result. 

It is interesting to examine whether one could obtain reliable estimates of the universal 
constants Ooo and p^o using the above window and to see whether such estimates can be 
tested against known results in the literature. It appears that the idea of such a test 
has not been tried in previous studies. In contrast, the main conclusion of the recent 



paper of Hilfer et al. 



2l| is that "the universal scaling function for the order-parameter 



distribution cannot be considered to be known from numerical simulations at present". 
Furthermore, it was stated in the conclusion of this paper that "even for the square Ising 
model a numerical estimation of the universal PDF requires sizes that are beyond present 
day computer recourses". Therefore, one would be tempted to think that, also the universal 
constants a^o and poo will approach their values slowly in the asymptotic limit. Thus, a 
reliability test for the values of Ooo and p^o estimated in the above a;-window will be useful 
even in the case where it fails. On the other hand, a test yielding a good comparison will be 
a strong verification of the proposal that the observed stable x-window can be used for the 
extraction of the asymptotic behavior of the tail regime of the universal scaling function. 

In order to construct a reliability test as discussed above, we shall now define following 
Bruce [l5[ the universal function Tiy) 



J^{y) = In 



(6) 



The limiting behavior of the function T{y) is controlled by the large- a; behavior of p{x) and 
as shown by Bruce [l5| one can find, by assuming the validity of the full ansatz ([5]), a large 
^/-expansion of J^{y) which reads as 



^(^/)^6^yi+V^ + iln 



oo 



(7) 



_ac^5{5 + 1) 

where h^o is a constant. The constant term in Eq. ([7]) is the universal Privman- Fisher 



coefficient [15|]. Its value for the 2D square Ising model is Uo = —0.639 912 [15|, l50|, and 



as shown in Ref. 



15| is related via the ansatz ([5]) to the universal constants Ooo and poo as 
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follows 



oo 



The expansion ([Tj) above has been fully illustrated by Bruce j.lSj], using appropriate y- 
windows to observe the development of the effective Uq (U^-^-^) in the large y-range. The 
agreement of this development towards the universal Privman-Fisher coefficient was shown 

n 

to be excellent (see Fig. 1 in Ref. |15l|). 

In Fig. we reproduce Fig. 1 of Ref. [l5|, using our numerical data for sizes L = 80 
and L = 120. The approach to the universal value Uo is again excellent and the L = 120 
data produce a slightly faster approach, as can be seen from this figure. The estimates of 
U^^^ in the range of — 5 deviate from the exact value by less than 0.2%. This is an 
explicit confirmation that a CrMES restriction is adequate for (numerical) studies of the 
tail regime. As pointed out by Bruce [isl , the numerical estimation of Uo via the numerical 
integration (see Eq. ([6])), illustrated above, is more reliable than attempting a direct Monte 
Carlo determination of free energies 1^, [o] or attempting to estimate Uo via the factors 
appearing in Eq. ([8]). In accordance with our earlier discussion, an accurate estimation of 
the universal constants a^o and poo would not be normally expected at these lattice sizes 
and corrections to scaling may complicate the situation. 

We proceed to describe an estimation procedure for the universal constants Ooo and poo, 
based on the stable window x = 1.2 — 1.3. As can be seen from the fitting attempts in Fig. HI 
the estimates for the critical exponent ip suffer from statistical fluctuations and are sensitive 
to both the DOS statistical method (WL or BH) and to the lattice size. It is therefore 
very important to repeat the statistical sampling and to use several lattice sizes in order to 
be able to observe systematic dependencies that may effect the estimated values. A fully 
converged implementation of the WL algorithm will be essential for the accurate application 
of the scheme and, of course, as usually repeated applications should improve the scheme. 

Thus, we have applied a well-saturated (see the discussion in Sec. [Tll) CRMES-WL(N- 
fold: 16-30) entropic scheme using a separation refinement 5 = 8 for the N-fold levels j = 
16 — 26 and a separation refinement S = IQ for the N-fold levels j = 27 — 30. In each 
run we used 4 such WL random walks to estimate the average WL DOS and to obtain 
the corresponding BH DOS. At the end of the run the cumulative Hwl{E, M) histograms 
were used with the above DOS's to obtain the corresponding universal PDF's. This process 
was repeated 4 times for each of the lattices sizes L = 90, 100, 110, 120 and 140. For each 
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lattice size, and each of the above runs (consisting of the 4 WL random walks) we fitted 
the function fl5a|) in the stable window x = 1.2 — 1.3 by fixing the value of the exponent 
ip to the expected value ip = 7. Thus, for each lattice size we obtained mean values and 
statistical errors for the following parameters: the universal constant a^o and poo, where Poo 
absorbs the identification of all peak-heights to 1. The multiplying factor, thus absorbed in 
our notation by using poo, is the peak-height p{x*) {poo = poo ■ pi^*)). The peak-height was 
also estimated for each run separately as well as the parameter A = L~^/®/ a/ (m^) which 
connects the universal PDF in the scaling variable x = m/ \/ {m?) jisjl with the universal 
PDF in the scaling variable z = rriL^^^ used by Hilfer et al. 19| {z = X~^x). Our fitting 
attempts were also repeated in this ^-representation using the equivalent stable window 
which has been taken to be z = 1.26 — 1.36. For the parameters in this 2;-representation the 
following notation was used: Pz,oo and a^^oo for the universal constants controlling the tail 
regime, and Pz{x*) for the peak-height. 

Using the above described fitting attempts we observed no systematic L-dependency, in 
the range L = 90 — 140. This observation concerns all the estimated parameters, and it 
appears that not even weak L-corrections should be applied. In fact the picture, obtained for 
each parameter, resembles effectively a statistical fluctuation around a mean value. There- 
fore it seems sensible to average the values of all parameters over the range L = 90 — 140 
(assuming no systematic L-dependency) and take as respective errors three standard devi- 
ations of the averaging process. The values of all parameters obtained from this hypothesis 
are collected in Table [T] and are presented for the methods corresponding to: the WL DOS 
(WL) and the BH DOS (BH). Both x- and ^-representations were used in the fittings and the 
average for A is also given in the footnote of Table [B Using these values in conjunction with 
Eq. ([8]) and the exact value for the state exponent 6 = 15 (note also that poo = Poo ■ p{x*)), 
the respective estimate for the universal Privman-Fisher coefficient is determined for each 
case and is presented in Table [B 

Discussing Table [T] we first observe that all parameters involved in the estimation of the 
universal Privman-Fisher coefficient have reasonable relative errors of the order of 2 — 4%. 
The corresponding relative errors for Uo are of the order of 4 — 10%. However, the values 
obtained for Uo are lower than the exact value and the underestimation is of the order of 7%. 
This, compared to the excellent agreement obtained by the numerical integration method 
shown in Fig. O and discussed above reveals the superiority of the integration method. The 
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origin of the present underestimation is not clear. One may think that, moving to a x- 
window corresponding to larger values of x could improve the above picture. Nevertheless, 
we have checked that this is not so, and that statistical errors become dominant as we move 
to larger x-values, making the estimation scheme unreliable. Finally, let us compare the 
constant a^^oo = 0.026... of Table [I] with the coefficient of the exponential term given in 
Eq. (7) of Ref. il9i]. This latter coefficient has the value K^l)^^^ = 0.023738... for the 
square Ising model. These two values are of the same order, although corresponding to 
different scaling limits 19j, and this is an indication in favor of the reliability of the stable 
x-window discussed here. It also strongly supports the proposal of Hilfer and Wilding [l9:i] 



to calculate critical finite-size scaling 
phase transitions developed by Hilfer 



unctions via the generalized classification theory of 



29 



30, 



31 



32|. It will be interesting to uncover the 



reasons behind the minor underestimation of the universal Privman-Fisher coefficient by the 
direct estimation of the universal constants a^o and p^o) as attempted in this work. 

Finally, we briefly discuss the behavior of the left tails of the critical distributions. This 
behavior has been considered by Hilfer et al. 2l|] by using rescaled probability densities 
Po{x)-, where x and Po are deflned in such a way that yield mean zero, unit norm, and 



unit variance 



2l[ |. This new scaled variable x should not be confused with the scaled 



variable x deflned in Eq. and appearing in Figs. [T] - HI These new distributions are 
appropriately translated with respect to the position of the peak of the distribution {xpeak) 
for each tail separately. For instance for the left tail one deflnes: Poi{x) = Po{xpeak — x) with 
X < Xpeak- Furthermore, these authors used plots of the functions q{y) = '^^°^^7og'x^°'^ versus 
y = logx, and compared their tail behavior with the standard Gaussian behavior. The large 
x-range of these plots is convenient for illustrating the possible developments of fat stretched 
exponential tails of the form ~ exp[— a;'^]. Using their method we have also tried to clarify the 
behavior of the left tails of our critical distributions for all sizes up to L = 140. Qualitatively, 
we found the same behavior with that of Ref. 2l|. In particular a data collapse for all sizes 
up to L = 140 was observed and the picture is very similar with their behavior (see: the right 
upper row of Fig. 6 in Ref. 2l|). However, it should be noted that, the data corresponding 
to the left tails suffer from relatively stronger fluctuations when compared with the data 
of the right tails. This is due to the inevitable more limited magnetic sampling in the left 
tails. Note that, the left tails are strongly influenced from a part of the energy spectrum in 
which a much larger number of spin conflgurations exists, compared with the lower energy 



15 



TABLE I: Estimates for the universal parameters and the Privman-Fisher coefficient. 



Method" Uo 





Poo 


p{x*) 


Poo 






WL 


0.6258(240) 


1.309(51) 


0.8190(449) 


0.0571(32) 


-0.589(63) 


BH 


0.6172(110) 


1.303(30) 


0.8042(234) 


0.0561(20) 


-0.599(34) 




Pz,oo 


Pz{x*) 




dz,oo 




WL 


0.4476(269) 


1.252(52) 


0.5604(156) 


0.0266(9) 


-0.589(33) 


BH 


0.4462(132) 


1.246(30) 


0.5560(213) 


0.0268(5) 


-0.599(40) 



WL: A = 0.951(11), BH: A = 0.953(5). 



part of the spectrum determining the right tails. However, using relatively larger x-windows 
(corresponding to at least 250 different magnetization values for L = 120) we have tried 
to repeat the small x-windows fitting attempts, as those shown in our Fig. [21 In this case 
we used as a test function the simple exponential behavior Poi{x) = qexp[—bx'^], and the 
total fitting range of the new scaled variable x was x = 1 — 3.25, which we assume describes 
the left tail regime. Note that this range corresponds to the range x = 0.87 — 0.28 of our 
Fig. [2] for the scaled variable defined in Eq. (jl]). For the lattice L = 120 the total fitting 
range includes about 2400 different magnetization values. Remarkably, we have obtained 
almost the same behavior for L = 50 and L = 120. For both sizes we found that for the 
above rather large range (1 — 3.25) the exponent uj steadily decreases (almost linearly for 
L = 120) from the value u = 0.75 to the value 0.15 passing through the value u = 0.5 in 
the neighborhood of x = 1.75. For larger values of x outside the above range (1 — 3.25) the 
statistical fluctuations do not permit reasonable fits. The described behavior is consistent 



with that of Ref. 



21| and since it is observed for the small and the larger lattice sizes we 



may suppose that it represents also the true asymptotic behavior. As a final remark we 
point out that the fat stretched exponential left tail ~ exp[— reported in Ref. 2l|] at 
the low temperature T = 1.5 can not be studied by our data. Such a study would require 
the application of the present scheme in a lower part of the energy space. 
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IV. CONCLUSIONS 



In the present Monte Carlo study we applied the CrMES-WL entropic samphng scheme, 
based on the high-levels of the Wang-Landau process in dominant energy subspaces, in 
order to obtain the order-parameter universal PDF of the square Ising model for large 
lattices. The efficiency of this method enabled us to clarify the asymptotic tail behavior 
of the universal distribution and to obtain reliable data for the universal parameters. In 
particular, we found that there exists a large stable window of the scaled order parameter in 
which the full ansatz for the tail regime is well obeyed. In a second stage, this window was 
used to estimate the equation of state exponent S and also to observe the behavior of the 
universal constants implicit in the functional form of the universal PDF and to approximate 
for the first time their values. The estimates of the universal constants appear to be reliable 
to within 3 — 7% of statistical errors. The excellent accuracy obtained for the universal 
Privman-Fisher coefficient, by appropriate numerical integration, was also illustrated and 
consist a concrete reliability test of the accuracy of our numerical scheme. 
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FIG. 1: Universal PDF's for a lattice of linear size L = 140, with a logarithmic scale on the vertical 
axis. Three different subspaces i?i,i?2, and R3 are used to obtain the corresponding curves. The 
wider subspace Ri {ie = 1950—3800) yields the PDF shown by the solid line. The dots demonstrate 
the PDF corresponding to the subspace R2 {ie = 2228 — 3514) and as can be seen from the inset 
the two cases coincide in their common part. R2 is the CrMES defined by the specific heat's 
accuracy condition, as discussed in the text. Finally, the dashed curve, i?3-case {ie = 2620 — 3800), 
illustrates that errors will be introduced by an inappropriate restriction of the energy space. 
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FIG. 2: Illustration of the universal scaling function p{x) for L = 60 and L = 120. The inset is an 
enlargement of the right tails for the case L = 120 and L = 140. Note that the peak-heights have 
been set equal to 1. In both cases a logarithmic scale on the vertical axis has been taken. 
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FIG. 3: Behavior of estimates of the exponent if)/! for lattice sizes L = 80, 120 and 140, for both 
the WL and BH methods used. The window x = 1.2 — 1.3 appears to be very stable for the sizes 
L = 120 and L = 140, where the estimates for the exponent ip remain close to the value ip = 7. 
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FIG. 4: Illustration of the three-paxameter fitting attempts for L = 100 and L = 120 using both 
WL and BH approximations of the universal PDF. The values of the estimates for the exponent 
and their fitting errors manifest the accuracy of our results and also the stability of the fitting 
attempts at the window x = 1.2 — 1.3. The small differences of the exponent estimates from the 
two different DOS's, within the same entropic sampling runs, reflect the sensitivity of the fitting 
attempts to statistical errors, growing with the lattice size. 
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FIG. 5: Variation of Uo^^ ■, for lattices L = 80 and L = 120, obtained as intercepts of linear fittings 



of Eq. ([7]) in small windows in the variable The y-m values correspond to the center of the 

fitting windows. The inset shows the linear fit according to Eq. ([7]) in the large y-part {L = 120). 
The dashed line marks the exact value of the amplitude Uq for the 2D Ising model i5Q]. approached 
by the corresponding intercepts. This figure should be compared to Fig. 1 of Ref. 
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